Loop integration results using numerical extrapolation for a non-scalar integral 

E. de Doncker a *, Y. Shimizu b , J. Fujimoto b , F. Yuasa b , K. Kaugars a , L. Cucos a , and J. Van Voorst a 

a Department of Computer Science, Western Michigan University, 
Kalamazoo, MI 49008, U. S. A. 

b High Energy Accelerator Research Organization (KEK), 
Oho 1-1, Tsukuba, Ibaraki, 305, Japan 

Loop integration results have been obtained using numerical integration and extrapolation. An extrapolation to the limit 
is performed with respect to a parameter in the integrand which tends to zero. Results are given for a non-scalar four-point 
diagram. Extensions to accommodate loop integration by existing integration packages are also discussed. These include: using 
previously generated partitions of the domain and roundoff error guards. 
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1. INTRODUCTION AND BACKGROUND 

In a general form, loop integrals used for cross sec- 
tion corrections are given by 
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where N is the number of propagators, L the number 
of loops, the momentum on the ^-th internal line is kg 
and the corresponding mass is mt.,1 < I < N. 

As a special case, scalar one-loop integrals of the 
form (-l)'7(167r 2 )/ n where 

In = /<S„_i (D„(x)-ie)"-^ dx ( 2 ) 

are obtained from Q by introducing Feynman param- 
eters and integrating over the loop momentum I. The 
integration region S n -i is the n— 1 dimensional unit 
simplex. 

For the simplest cases, the results can be obtained 
analytically. Numerical techniques have been suc- 
cessful with considerable analytic manipulation (see, 
e.g., 1 1121 ). In previous work Q, we reported results 
for integrals of the form treated numerically us- 
ing extrapolation by the e-algorithm |4 1. We will now 
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consider the case of a one-loop integral where the nu- 
merator in the integrand is a polynomial of the Feyn- 
man parameters. A sample problem involving the 
e~ e + — > W~ W + interaction is given in the next sec- 
tion of this paper. Results for this problem are given 
in Section|5] Section|4]discusses enhancements to the 
Parlnt parallel integration package. 

2. NON-SCALAR INTEGRAL 

The matrix element of one-loop corrections is 
given by the real part of the product of a one-loop am- 
plitude and the (conjugate of) a tree amplitude. Fig- 
ure |2 shows an example of a box diagram and a tree 
diagram of a Z-boson exchange for the interaction 
e~e + — > W~ W + . The Feynman diagram and the 
corresponding matrix element are generated automat- 
ically by GRACE-loopO system. 

After introducing the Feynman parameters as in 
Figure |2 and integrating over the loop momentum, 
the matrix element is of the following form, 
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Figure|2]shows the D 4 = surface of the singular- 
ity over —1 < x, y < 1, and delineates the integration 
domain S3 . 
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Figure 1 . Feynman diagram for e e + — > W W + 



f and g are polynomials of Feynman parameters, 
of which the coefficients are determined by exter- 
nal momenta and masses of internal lines. Here 
M z = 91.187GeV,Miv = 80.22GeV,m e = 
0.511MeV, y/s = 500 GeV and 9 = 
^(Pe-iPff-)- Th e numerical results are evaluated 
for cos 9 = 0.956811390. 

The generalized non-linear gauges|5) are imple- 
mented for the amplitude. The result depends on 
the gauge parameters because only one diagram is 
picked up. For the numerical evaluation, the non- 
linear gauge parameters are set as a = 2, $ = 3, 5 = 
4, e = 5 and k = 6. 

3. GRAPH 216 RESULTS 

Table [2 illustrates the use of the e-algorithm for 
the integral computation of the term involving / (the 
symbolic code of which has about 2000 lines as FOR- 
TRAN code). We show the results of the extrapo- 
lation for the real part of Mi(f, 0; e). The method is 
based on generating a sequence of integral values cor- 
responding to a geometric sequence of e and extrapo- 
lating to the limit as e — > 0. 

The table shows the sequence of integral approx- 
imations for e = 1.2 30- '', I = 0, 1, . . . (obtained 
numerically) in the first (leftmost) column. Using the 
integral approximations corresponding to £ = 0, 1, 2, 
the first extrapolated result is obtained (top element 
of column 2). Using the £ = 3 element of column 




Figure 2. D4 = surface 



1, the second extrapolated result is obtained in col- 
umn 2; the £ = 4 element of column 1 is then used 
to generate the third element of column 2 and the top 
element of column 3. In all iterations following, the 
new element in column 1 is used to generate a new 
lower diagonal of the triangular table. 

The table elements are shown to 8 -digit accuracy, 
which is the final accuracy obtained in this run. Con- 
vergence is apparent down the columns and along the 
lower diagonal. Relying on a heuristic error estimate 
of the table elements along the lower diagonal, an el- 
ement is selected as the result (printed boldface). The 
result calculated analytically is -0.647837287. 

To generate the integral approximation in the first 
column, we used an iterated integration where the 
adaptive Quadpack |6) routine DQAGE was used 
in each direction, requesting a relative accuracy of 
10~ 10 . So far, this technique has outperformed other 
numerical integration approaches using multivariate 
(cubature) rules. 

4. PARENT ENHANCEMENTS 

Parlnt is a software package for parallel multivari- 
ate integration f2J. It has components for multivari- 
ate integration using Monte Carlo (MC), Quasi-Monte 
Carlo (QMC) and adaptive methods. Parlnt is written 
in C and runs over MPI 1 8 1 on a distributed platform. 
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Figure 3. 2D of ID integrand for Graph 216 real part 



4.1. Iterated integration 

While the adaptive approach could not be applied 
directly using 3D multivariate rules, results to 6-figure 
accuracy were obtained by treating the problem as a 
2D integration of a ID integral. The ID inner integral 
was calculated with Quadpack routine DQAGE. 

The 2D integration was performed with Parlnt and 
with its Fortran sequential predecessor, DCUHRE (9). 
The local region error estimate was changed to make 
it less conservative. Figure [3] illustrates the integrand 
of the 2D problem for e = 1.2 25 , which was drawn 
using evaluation points of the integration. We are cur- 
rently considering a design of Parlnt which will al- 
low incorporating iterated integration in a transparent 
way. 

It recently came to our attention that in the work by 
Binoth, Heinrich and Kauer iTOI . 3D box integrals are 
obtained by performing the inner integration analyti- 
cally, which leaves the resulting 2D integrand with an 
integrable (though still problematic) singularity. Note 
further that their 3D box together with 2D vertex di- 
agram evaluations are at the basis of reductions per- 
formed to treat scalar hexagon integrals. 

4.2. Re-use of subregions between extrapolations 

A sequence of extrapolation steps uses a series of 
similar integrations which share similar subregions 
when performed adaptively. At each step of the ex- 
trapolation, ParInt can avoid a significant number 
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of region evaluations by re-using previously "discov- 
ered" subregions as the initial set of regions for the 
next extrapolation step, potentially avoiding a signifi- 
cant amount of computation. 

ParInt has been modified to support this activity 
by storing the active integration regions at the end of 
every parallel run and providing the user with the op- 
tion to load a set of regions to initialize a run. Re- 
gions may be saved either locally on each compute 
node or in a single file managed by the controller. Re- 
gions loaded at the start of a subsequent computation 
may also be read from a single global file or individual 
files on each compute node. We are currently devel- 
oping a distributed I/O system which will allow com- 
pute nodes to retrieve previously saved regions from 
files on any other compute node II II . 



4.3. Kahan summation 

The global adaptive integration algorithm first de- 
veloped by De Ridder and Van Dooren 1121 is also 
used by ParInt. At each step, one region (per worker) 
is selected and subdivided into subregions. The se- 
lected integration rule is applied over each subregion. 
Next, the estimated error and result for the selected 
region and subregions must be subtracted from and 
added to the total estimated error and result, respec- 
tively. For difficult problems, ParInt will select many 
regions and subdivide them. Numerical summation 
of millions of terms can introduce round off error and 
greatly reduce the accuracy of the result and estimated 
error in a numerical integration routine. 

We have looked at several techniques to reduce 
round off error in sums with a large number of terms. 
Each of these techniques has its own merits and flaws. 
A good method would be one whose accuracy does 
not depend on the number of terms in the sum and 
would not greatly impact the runtime performance of 
a numerical integration routine. A compensated sum- 
mation method developed by W. Kahan 111 31 and fur- 
ther studied by N. Higham 1141 best fits these needs. 
Several advantages of this method are low computa- 
tional overhead, low storage requirements, and in er- 
ror analysis it is shown to have an error constant of 
order 1. 



5. CONCLUSIONS 

We presented results for a non-scalar one-loop box 
diagram, where the integral is obtained using numeri- 
cal integration and extrapolation with the e-algorithm. 
We described enhancements to the ParInt parallel in- 
tegration package, which are in various stages of de- 
velopment. Furthermore, in future work, we plan to 
investigate combinations of our numerical methods 
with symbolic techniques. 
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